function [x,normk,it,nrm] = solverNT_mc(in,xo)
% this function solves the coarse mesh problem via what is essentially
% a classical newton method.  Here, we can compute nearly completely the
% Jacobian given the definition of the nonlinear function.  Only values
% which are explict functions of keff need to be approximated via finite
% differences.  For reference, the nonlinear function is:
%  FF = [ R*M-x(end)*eye(length(R(:,1)))  zeros(length(R(:,1)),1) 
%         F*M                            -sum(A*M*x(1:end-2))-L          ];  
% z1 = FF*(x(1:end-1)); 
% z2 = 0.5 - 0.5*sum(x(1:end-2));

disp('-------------------------------------------')
disp('--      beginning newton iterations      --')
disp('-------------------------------------------')

normk = 0;
x           = xo;%;[xo' 1.0]'; % add in current eigenvalue
x(end)=1;
%x(end-1)=0.98988;
xx(:,1)     = x;
[z,R,F,A,L] = resid2D(x,in); % initial residual
itmx        = 20; 
it          = 1;
nrm(1)      = norm(z);

while ( nrm(it) > in.tol(1)*nrm(1)+in.tol(2) && it < itmx )
    disp(['  iteration number: ', num2str(it), ' norm(x) = ', num2str(nrm(it))])
   % fp          = jacob2D(x,in,z); % FD appx (for verification)
    fp          = jacob2Da(x,in,z,R,F,A,L); % hand coded
    s           = -fp\z;
    x           = x + s;
    it          = it+1; 
    xx(:,it)    = x; 
    [z,R,F,A,L] = resid2D(x,in);
    nrm(it)     = norm(z);
end

end % solverNT

function [z,R,F,A,L] = resid2D(x,in)
%% function z = resid2D(x,input)
% This function computes the nonlinear residual
% Here, x(1:end-2) corresponds to the legendre coefficients
%       x(end-1) is keff
%       x(end)   is the current eigenvalue (lambda)
% NOTE: I return the most recent response functions
M     = in.M;
order = in.order;
numel = in.numel;
numg  = in.numg;

[R,F,A,L] = combmesh_center_2d_resp(in,x);
loss = A*M*x(1:end-2)+L;%/x(end);

%compute the nonlinear function, FF(x)*x = 0
FF = [ R*M-x(end)*eye(length(R(:,1)))  zeros(length(R(:,1)),1)
       F*M                            -loss  ];
z1 = FF*(x(1:end-1));
z2 = 0.5 - 0.5*x(1:end-2)'*x(1:end-2);
z  = [z1' z2']';

end % resid2D

function fp = jacob2Da(x,in,z,R,F,A,L)
%  x  = current unknowns vector
%  in = problem data
%  z  = last computed residual
M = in.M;
n   = length(x);
eps = 1.e-8;
fp  = zeros(n);

%[R,F,A,L] = combmesh_center_2d_resp(in,x);

fp(1:end-2,1:end-2) = R*M - x(end)*eye(n-2);
fp(n-1,1:n-2)       = (F-x(end-1)*(A))*M-x(end-1)*L;
fp(n,1:n-2)         = -x(1:end-2)';
fp(1:n-2,n)         = -x(1:end-2);

% compute the perturbed response functions
x(end-1) = x(end-1)+eps;

[R2,F2,A2,L2] = combmesh_center_2d_resp(in,x);

Rk = (R2-R)/eps;
Fk = (F2-F)/eps;
Ak = (A2-A)/eps;
Lk = (L2-L)/eps;
fp(n-1,n-1)     = (Fk-A-Ak*x(end-1))*M*x(1:end-2) - (L + Lk)*x(end-1);
fp(1:n-2,n-1)   = Rk*M*x(1:end-2);
end % jacob2D

function fp = jacob2D(x,in,f0)
% fp = jacob(f,x,fo)
%    this function produces a crude forward finite difference
%    approximation for the function f(x) about x; fo is the initial
%    evaluation of f
n   = length(x);
%f0  = feval(f,x);
eps = 1e-8;
fp  = zeros(n);
for j=1:n
    v           = zeros(n,1); 
    v(j,1)      = 1.0;
    xx          = x+v*eps;
    f1          = resid2D(xx,in);
    fp(:,j)     = (f1 - f0)/eps;
end
end

function [R,F,A,L] = combmesh_center_2d_resp(in,x)
M = in.M;
order = in.order;
numel = in.numel;
numg = in.numg;
% initialize response matrices
R = zeros( numel*numg*(order+1)*4 );
F = zeros( 1, numel*numg*(order+1)*4 );
A = zeros( 1, numel*numg*(order+1)*4 );
L = zeros( 4, numel*numg*(order+1)*4 );
lenR = numg*(order+1)*4;
% build mm
for i = 1:length(M(:,1))
    if ( sum(M(:,i)) + sum(M(i,:)) )==0
        mm(i)=1;
    else
        mm(i)=0;
    end
end
mm = diag(mm(1:numg*(order+1):end));
%% Get response matrices
[Rb Fb Ab Lb]  = mesh_center_2d_resp(x(end-1),in);
% this gets the info for each element type; now, need to place
I = length(in.elements(:,1));
J = length(in.elements(1,:));
k = 0;
for j = 1:J      % rows
    for i = 1:I  % cols
        if in.elements(i,j) > 0
            % I am the kth element
            k = k + 1;
            type = in.elements(i,j);
        else
            % else I'm just void
            continue
        end
        R( ((k-1)*lenR+1):(k*lenR),((k-1)*lenR+1):(k*lenR) ) = Rb(type,:,:);
        F(1, ((k-1)*lenR+1):(k*lenR) ) = Fb(type,:);
        A(1, ((k-1)*lenR+1):(k*lenR) ) = Ab(type,:);
    end
end
R = sparse(R);
tmp = numg*(order+1);
k=0;
%% Get the leakage
for j = 1:J      % rows
    for i = 1:I  % cols
        if in.elements(i,j) > 0
            % I am the kth element
            k = k + 1;
            type = in.elements(i,j);
        else
            % else I'm just void
            continue
        end
        
        myL(:,:) = Lb(type,:,:);
        my1 = myL(1,:);
        my2 = myL(2,:);
        my3 = myL(3,:);
        my4 = myL(4,:);
        % for an element, we have
        Ltmp = [ my1 my2 my3 my4
                 my2 my1 my4 my3
                 my4 my3 my1 my2
                 my3 my4 my2 my1 ];
        % build the L matrix
        LK( (k-1)*4+1:(k-1)*4+4, (k-1)*4*tmp+1:k*4*tmp ) = Ltmp;
    end
end
% NOTE: 3.16.2010 -- I added the M because the input is actually J+
%  and LK is meant to couple with J-
L = sum( mm*( (M*x(1:end-2))'*LK')' );
end % combmesh_center_2d_resp
